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Abstract 



We present a method for sampling microscopic configurations of a pliysical system dis- 
Q\ • tributed according to a canonical (Boltzmann-Gibbs) measure, with a constraint holding in 

e average. Assuming that the constraint can be controlled by the volume and/or the tempera- 

ture of the system, and considering the control parameter as a dynamical variable, a sampling 
strategy based on a nonlinear stochastic process is proposed. Convergence results for this dy- 
namics are proved using entropy estimates. As an application, we consider the computation 
^ . of points along the Hugoniot curve, which are equilibrium states obtained after equilibration 

of a material heated and compressed by a shock wave. 



Statistical physics provides a way to obtain macroscopic quantities starting from systems de- 
scribed at the microscopic level. In this framework, the state of the system is described by some 
probabiUty measure, the precise choice of the measure depending on the invariant quantities at 
Q ' hand. A fundamental statistical ensemble for instance is the microcanonical ensemble which is 

O . generated by the ergodic Hmit of the Hamiltonian dynamics, and corresponds to isolated systems. 

Another classical ensemble is the canonical ensemble, for which the temperature is constant. In 

fT^ , one possible derivation of the associated probability measure, the statistical entropy is maximized 

^ ' subject to the constraint that the average energy of the system is fixed, and the temperature is 

OO , related to the Lagrange multipHer of the former constraint [3] . More generally, it may be of interest 

to consider ensembles such that some constraint (not only the energy) is satisfied on average. The 
challenge is to have an expUcit description of the corresponding thermodynamic ensemble in terms 
of the associated probability measure. In an effort to define the statistical distributions sampled 
by methods aiming at satisfying constraints on average, we first reformulate here the problem of 

QQ ' sampling with constraints satisfied in average in the context of the canonical ensemble. 

^^ , The paper is organized as follows. In Section[Tl we formulate precisely the problem of sampling 

canonical ensembles with a constraint fixed in average, giving a specific example and recalling 
some previous sampling strategies. In Section [51 we present a nonlinear stochastic process with a 
feedback term depending on the expectation of the constraint, who admits the target measure as 
a stationary state. We also prove the convergence to this target measure using entropy estimates 
for initial data not too far from the limiting state (the proofs being presented in Section [3|) . In 
Section [3l we discuss the numerical implementation, proposing a single trajectory discretization, 
and use the method to compute the Hugoniot curve of Argon. 

1 Mathematical formulation of the problem 

Consider a microscopic system of N particles in a space of dimension d = 3, described by their 
positions q = (qi, . . . .qn) and momenta p = (pi, . . . ,pn), with associated mass matrix M = 
Diag(77ii, . . . , ttt-at), and interacting through a potential V. Periodic boundary conditions are used, 



and the particles then stay in a domain T) = L^T x LyT x L^T (T denoting the one-dimensional 
torus R/Z). The volume of the domain is denoted by \'D\. The phase-space Q, is 

The central quantity describing the system is the Hamiltonian 

N 2 
i—l 

The canonical measure associated with the Hamiltonian ([T]) has a density [1] 

7rp.T(g,p) = ^c-^"^'i'P\ /3-1 = fceT, (2) 

Zix>.T 

where fee is the Boltzmann constant, and the partition function Zxi,t is a normalization factor so 
that ^ is indeed a probability measure: 

Zt,.t^ I c-^^^-^'PUqcip. (3) 

We have indicated explicitly the dependence of the canonical measure ^ and the partition func- 
tion Q on the temperature T and the domain V. Average thermodynamic properties of the 
system can be computed as averages of functions of the microscopic variables (the so-called ob- 
servables) with respect to the canonical measure at a temperature T and for a given simulation 
boxP: 

{A)v.T = / A{q,p)-Kv,T{q,p)dqdp. (4) 

Canonical simulations are performed at a fixed volume. However, this is an ill-defined notion since 
there are infinitely many simulations domains having the same volume. In practice, the shape of 
the simulation box (for instance cubic) is kept constant, and the volume is changed using a scaling 
of the positions q. We will in any case assume in the sequel that the knowledge of the volume |I?| 
is enough to characterize the domain T). 

1.1 Equilibrium sampling 

In most cases, equilibrium sampling at a fixed temperature and for a fixed geometry T) is performed 
(see in [HIT] for references on sampling methods), and the output of the simulation is the average 
property {A)xi,t- This can be done for instance by means of ergodic results for the Langevin 
dynamics 

dqt = M^^ptdt, 

dpt = -WViqt) dt - £,M-^pt dt -f ^2£,k^TdWu ^^^ 

where ^ > is physically interpreted as a friction, and Wt is a standard 3A^-dimensional Brownian 
motion. Under mild assumptions on the potential V [7], 



lim -/ A{qt,pt)dt = {A)t:.t a.s. 



for almost all initial conditions. When the observable A considered depends only on the positions 
of the particles (which is not too restrictive since the momentum contribution of many observables 
can be explicitly computed), the overdamped Langevin dynamics 



dqt = -VViqt) dt + ^2knTdWu 
may be used as well since (again, under some assumptions on the potential [7]) 
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lim - / A{qt)dt = '-^ = {A)v,T a.s. 



1.2 Satisfying constraints in average 

The inverse sampling problem is sometimes of interest: given some value Aq of an average property 
(for instance, some average energy), what should the temperature and/or the volume of the system 
be? Answering this question allows to sample configurations of the system canonically distributed 
and satisfying the constraint A{q,p) ~ A^ in average (in the sense that {A)x>,t — Aq). In the 
sequel we will focus on constraints depending on T, and therefore drop the mention of the volume 
in the notation of the canonical averages when it is not relevant. Constraints on the volume can 
be handled similarly. Moreover, upon replacing Ahy A — Aq, it can be assumed that Aq = 0, so 
that the problem under investigation is 



Find T such that {A)t = 0. 



(6) 



In order to have a well-defined problem (and possibly upon replacing A by —A), it is assumed 
that 



ure 



Assumption 1. There exists an interval I^ = \T^-^^,T^^^ (with TA^,, > 0), a temperat 
T* £ (T^i^^T^^^), and constants a,a > such that 

and 

yrj.^ rA ^ (^)t - {A)t* . 

yr e It, a< y_y^ < a. 

The assumption on the observable is satisfied as soon as T f-> {A)t is locally C^ around a value 
T* such that {A)t' = and dr {{A)T)\rp, > 0. An example is illustrated in Figurc[2]for the test 
case considered in Section [ 

1.3 Some examples 

1.3.1 Hugoniot curves 

The method presented in this paper is illustrated on Hugoniot sampling (see Section 13. 2p . The 
variations of macroscopic quantities across a shock interface are governed by the Rankine-Hugoniot 
relations, which relate the jumps of the quantities under investigation (pressure, density, velocities) 
to the velocity of the shock front. The third Rankine-Hugoniot conservation law for the Euler 
equation governing the hydrodynamic evolution of the fluid reads (macroscopic quantities are 
denoted by curly letters) 

£-£o-l{V + Vo)iVo-V) = 0. (7) 

In this expression, £ is the internal energy of the fluid, V its pressure, and V its volume. The 
subscript refers to the initial state, the other quantities are evaluated at a shocked state, after 
equilibration. The Hugoniot curve corresponds to all the possible states satisfying ([7]). In practice, 
it is computed by considering shocks of different strengths, inducing various compressions. 

A reference temperature Tq and a simulation cell 'Dq, for instance 2?o = (Lo^)^! characterize 
the equilibrium state before the shock. In numerical experiments, the compression rate 

\V\ 
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is varied from 1 to some maximal compression rate < Cmax < 1, so that T) — {c^''^LoT)^ 
when the compression is isotropic. Since all macroscopic quantities arising in the hydrodynamic 
equations are obtained relying on some local thermodynamic equilibrium assumption, ([7]) can be 
reformulated at the microscopic level using statistical mechanics. For a given compression rate c, 

(iJ)|^|,T - (i^)lPol.To - 1{{P)\V\,T + (0|Po|,To)(l^o| - m = 0, 



where the pressure observable for a domain T) is 

1 ^ 2 

Notice that we indicate exphcitly the dependence on the volumes for clarity, though the initial 
and final volumes are given. The final temperature T is the only unknown, and it satisfies 

{H{q,p) - (if)|p„|,T„ + \{P{q,p) + {P)\v„lT,){l - c)|Po|\ = 0. 

\ ^ / \'D\,T 

Introducing the Hugoniot observable (parametrized by the compression parameter c) 

A,{q,p) = H{q,p) - (if>|p„i,To + \{P{q,p) + (P> |P„|,T„)(1 - c)|2?o|, (8) 

the Hugoniot problem can be reformulated as 

For a given compression c,nax < c < I, find T such that {A^) c\'d^\t — 0- (9) 

The compression rate c parametrizes a curve in the (P, T) diagram, called the Hugoniot curve. 
Note that this curve does not correspond to a thermodynamic path. 

Since shock waves propagate in one direction (for instance parallel to the x axis) , anisotropic 
version of the Hugoniot problem are of interest. In this case, the compression acts in the x direction 
only, and V = cLqT x (LqIT)^. The average pressure P is replaced by the P^x component of the 
pressure tensor: 

where qi^xiPi,x denote respectively the x components of the position and momenta of the i-th 
particle, and the observable Ac is replaced by 

Axx,Mp) = H{q,p) - {H)iVo\.T„ + ^(Pxxiq^p) + {Pxx)\Vo\,ToK^ ~ c)\Vo\. 

The XX component of the initial pressure tensor {Pxx)\t>o\.To ma-y be replaced by the average 
pressure {P)\Vo\.To when the initial state is isotropic. 

1.3.2 A more general case 

More generally, Hamiltonians depending on some external parameter A can be considered, with 
associated obscrvablcs A controlled by the value of A. For instance, A could be the intensity of a 
magnetic field for a spin system, and A the total magnetization. The problem is then to find A 
such that 

/ A{q,p)e-l^"^^'i^P'^dqdp 
{A)t>.t.x = ^-^ = 0, 

/ e-^^^^^'PUqdp 
Jn 

the volume V and the temperature T being fixed. 

1.4 Previous sampling algorithms 

1.4.1 Newton's method 

An obvious sampling strategy to compute the zero of the function T i-^ F{T) ~ {A)t is to resort to 
Newton's algorithm. It is expected that such an approach converges in a few steps if the function 
F is well behaved. However, the derivative F'{T) often involves the computation of covariances, 



which are known to require long samphng times to be rehably computed. For instance, when 
trying to sample configurations with the average energy fixed, A ^ H and 

A way to bypass this difficulty is to estimate numerically the derivative by computing averages 
at two temperatures T + AT and T — AT around the temperature T, and to approximate the 
derivative as 

^ n A\ \ ^ (^>T+AT - (^)t-AT 
dT{{A)T) ^ ,^^ . 

This method is robust, and is indeed used in Monte-Carlo studies of the Hugoniot problem [S]. 
It requires however very carefully converged canonical samplings in order to have a numerical 
estimate of the derivative not polluted by sampling errors. 

1.4.2 New thermodynamic ensembles 

It would be better to have a more automatic procedure, consisting for instance of a single molecular 
dynamics trajectory. In this case, convergence needs only be checked once, and the simulation 
can be runned longer if needed. This was the motivation for the Hugoniostat method |15] . which 
can be extended to any constraint whose values are controlled by the temperature T. The idea 
of the method is to satisfy approximately the constraint A(q,p) = at all times. To this end, a 
convenient dynamics in the Nose-Hoover fashion [T71 [13] is postulated: 

q = M^^P, 

p = -VV{q)-ip, 
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where i/ > is homogeneous to a frequency, and A^^i is some reference value of the observable. 
When the instantaneous value of the observable is lower than the target value 0, ^ decreases 
and the friction is reduced, so that the (kinetic) temperature in the system, and therefore the 
average value of the observable, can increase; while the friction increases for too large values of 
the observable, and so, energy is removed from the system and the temperature decreases. This 
feedback process ensures indeed that values of the observable around the target value are sampled, 
though it is unclear how the temperature of the system is defined. For computational purposes, 
the temperature is defined as the average kinetic temperature along the trajectory. 

The numerical results obtained with the dynamics pUj) are correct for the Hugoniot problem 
since they are in agreement with numerical results from shock wave simulations, see [15j for more 
precisions. However, the question of the choice of the frequency v in (jlOp remains, as usual for 
Nose-like dynamics. More importantly, from a fundamental viewpoint, an unpleasant feature is 
that the thermodynamic ensemble {i.e. the probability measure on phase space) associated with 
the dynamics used is not known. It may be defined as the ergodic limit of the dynamics, for 
a given set of parameters and initial conditions. Even if this is possible, it is unclear whether 
the limiting measure depends on the parameters used and on the initial conditions, or not. The 
problem remains even if the Nose-Hoover part of the above dynamics is replaced by a Langevin 
dynamics (for which ergodicity results are usually easier to state) or any dynamics ergodic for the 
canonical measure (see for references on the convergence properties of some sampling methods) . 
For instance, the usual Langevin dynamics (O at temperature T may be extended to the case of 
constraints satisfied in average by considering the temperature as a dynamical variable: 



dqt = M ^ptdt, 



dqt = M pt at, 

dpt = -VViqt)dt-^M-^ptdt + VWk^tdWt, 



(11) 
dTt = -ly — Tretdt, 



A- 
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where i' > is again some frequency, and Tref, A^ci are reference values of the temperature 
and the observable respectively. When the instantaneous value of the observable is too low, the 
temperature is increased, while it is decreased when it is too large. For computational purposes, 
the temperature may be defined as the time average of Tt over a trajectory, though, as for the 
dynamics (jlOp , the meaning of this quantity is unclear since there is no obvious definition of the 
thermodynamic temperature within the ensemble generated by the dynamics pip . Therefore, pi|l 
is not satisfactory. 

2 A nonequilibrium sampling strategy 

2.1 General algorithm 

To solve the sampling problem ([6]), wc consider the temperature as a variable, so that the system is 
described by the variables {q,p, T). The idea is that the update of the current configuration {q,p) 
is governed by a dynamics consistent with canonical sampling at the instantaneous temperature 
T{t) (such as Nose- Hoover methods [T71[T3], Metropolis-Hastings algorithms [TBI ED, Langevin or 
overdamped Langevin dynamics,...), while the temperature is updated depending on the current 
estimate of {A)x(ty For instance, consider the following system when the underlying dynamics is 
of overdamped Langevin type: 



I dqt = -VV{qt)dt+^2kjiT{t)dWu 
\ T'{t) = -7lE(Afe)), 

dqt = M^ pt dt, 



dpt = -VViqt) dt - ^M-^pt dt + ^2^kBT{t) dWt, (13) 

T'(t) = -7E(A(g,,pO), 
when the underlying dynamics is of Langevin type. In both cases, 7 > determines the time 
scale of the temperature feedback term, and Wt is a standard SA'^-dimensional Brownian motion. 
To prove mathematical results, it is convenient to use (elliptic) overdamped Langevin dynamics, 
since the associated Fokker-Planck equation describing the evolution of the law of the process qt is 
parabolic. However, we will use the hypoelliptic Langevin dynamics in the numerical applications 
presented in Section [3] since this dynamics is usually a more efficient sampling device [7] . When 
the overdamped Langevin dynamics is used, the canonical measure ^t in position space has the 
density 

and averages with respect to the canonical measure jit arc still denoted {A)t = / A{q) ^xiq) dq. 

J M 

2.1.1 Motivation 

The dynamics (fT^ and p^ can be motivated as follows. If the configurations followed adiabati- 
cally the temperature changes, the positions qt (and possibly the momenta pt) would be distributed 
canonically at the temperature T{t) at all times, and 

r'(i) = -7(A)T(*). 

Assumption ([T]) then shows that T{t) — > T* . Indeed, starting from a temperature T(0) in the 
interval /^ of Assumption [l] it is easily seen that t i-^ {T{t) — T*)^ is a decreasing function: if 
T{t) G /^, then 



dt 



\{T{t)~T*f 



= -j{A)r^t){T{t) T*) = -7((A)t(*) - {A)T'){T{t) T*) 
< -7a(T(i)-r*)2 < 0, 



and so, T{t) £ I^ whenever r(0) G I^. Suppose for instance T(0) > T*. Then T{t) > T* at all 
times and, with Assumption [l] 

-Mnt) - T*) < T'{t) = -7((A)T(t) - {A)t') < ~ia{T{t) - T*). 

Gronwall's lemma then implies 

T* + (r(0) - T*) e""'^* < T{t) <T* + (T(0) - T*) e-"^*, 

which shows the convergence T{t) -^ T* as t ^ +00. 

Of course, the positions are not canonically distributed at all times. The hope is that, even 
if equilibrium is not maintained at all times, the dynamics can still converge if the typical time 
arising in the temperature update is small enough compared to the typical relaxation time of 
the dynamics, so that the spatial relaxation happens gradually, as T{t) approaches T*. This 
statement is made precise in Section [221 (see Theorems [2] and |3]) . From a practical viewpoint, 
the interest of such a strategy is that, if the dynamics is started off equilibrium (for instance 
at a temperature much larger than the target temperature T*), then it is not necessary to fully 
relax the system at the starting temperature (as is the case for instance for the Newton method 
presented in Section fl. 4. ip . and so, computational savings may be anticipated. The convergence 
proof is however only valid for initial data not too far from the stationary state. This is related to 
the fact that the temperature has to remain positive, and so, it is not possible in general to start 
arbitrarily far from equilibrium. 

Remark 1 (Possible extensions). A nonlinear temperature feedback may be considered as well: 

r'(t)=.-7/(E(A(gO)), (14) 

under some assumptions on the function f - see Remark\^for more precisions. 

It is also possible to use a time varying parameter ^[t) for the temperature update, for instance 
smaller at the beginning of the simulation, and then larger once the equilibrium state is approached. 
The proofs presented below may be extended to this case provided ^(t) remains in a properly chosen 
interval around an admissible value of "f (as given by Theorems\^ and\3\) . 

2.1.2 Partial differential equation reformulation 

The nonlinear partial differential equation (PDE) reformulation of the dynamics (fT^ is 



dti^ = fcBr(t)v 



tJ-Tit)"^ 



\l^T{t) 



(15) 



r{t) = -7/ A{q)^{t,q)dq, 

Jm 



where the periodic function q Cz A4 t-^ ip{t, q) is the law of the process qt at time t. 

The PDE (fT5|l clearly admits {T* , fiT") as a stationary solution. Since the first equation of 
(|15p is of parabolic type, standard techniques may be used to prove the (short time) existence 
of a unique solution, for initial conditions smooth enough, and under the following smoothness 
assumptions on the potential and the observable: 

Assumption 2. The observable A e C^(A^) and V £ C^(A^). 

In particular, the observable is bounded, which is important to ensure that the temperature 
remains positive. Indeed, if this was not the case, it would be possible to obtain negative tempera- 
tures arbitrarily fast by concentrating the initial density ip^ around singularities of the observable. 
The boundcdness of A ensures some delay in the feedback. 



Theorem 1. When Assumption\^ holds, and for a given initial condition [T^^ip'^), with T" > 
and 

V'" e H2(X), ^° > 0, f ip° = 1, 

JM 

there exists a time t > — — — — > such that (llSp has a unique solution (T, ip) with T G 
Ci([0,r],R) andil>£ C°([0, r],H2(A^)). Moreover, V > 0, and if: > when tp" > 0. 

The proof is presented in Section j4.1l The regularity assumptions on ip'^ is not too restrictive 
since an arbitrary initial condition ■0° G L^(A^) can be evolved using an overdamped Langevin 
dynamics at T° for some time ri„it > before turning on the temperature feedback term. Then, 
the resulting distribution ■0° = '(/'(rinit, ■) is smooth thanks to the regularizing properties of the 
parabolic Fokker-Planck equation. To obtain the uniqueness and existence of the solution at all 
times, it is necessary to show that the temperature remains isolated from 0. This is proved in 
Section [2?2l together with the convergence to the stationary state, under certain conditions on the 
initial entropy (see Theorems [5] and [3]) . 

2.2 Convergence results for initial conditions close to the fixed-point 

To study the convergence of the nonlinear dynamics (|15p under its PDE form, entropy methods 
can be used (see for instance the review papers [IHIE])- The total entropy of the system is the sum 
of a spatial entropy E{t), related to the distribution of configurations at time t, and a temperature 
term: 

£{t)^E{t) + \{T{t)-T*r, 
with 

Eit)= [ /l(/)/.T(*), ./= "^ 



M l^T{t) 

Notice that the spatial entropy E{t) measures the distance of the law of the process -0 to the 
instantaneous equilibrium measure fiT{t): ^nd not to some fixed reference measure. Two classical 
choices for the spatial entropy are the relative entropy distance (which corresponds to h(x) = 

xlnx — X + 1 > 0) and L^ distances {h(x) = —{x — 1)^). We present convergence results for both 

cases (see Theorems [2] and |31), which are very similar in spirit, but involve different functional 
settings for the initial conditions. In particular, the conditions on the initial condition are more 
stringent in Theorem [31 but the assumptions to be verified by the potential are less demanding. 

The convergence is ensured provided £{t) ^ as t -^ +c». This implies in particular T{t) -^ 
T* , and also 

||'0(t, O-Mt-IIli < ll'0(i, •) -AiT(t)llLi + llMT(t) - A^t-IIli < ChE{t) + \\^iT(t) -Mt-IIli ^ 0, 

where the constant Ch depends on the chosen entropy density h (using a Csizar-KuUback inequality 
for relative entropy estimates, and a Cauchy-Schwarz inequality for L^ entropies). The strategy 
of the proof is to obtain a Gronwall inequality for £(t). This ensures also that the temperature is 
well defined (in particular, T{t) > at all times). Indeed, if the total entropy is decreasing, then 
\T{t) — T*\ is bounded by ^J2£{Q), and the temperature remains isolated from for initial data 
whose total entropy is small enough. 

In the proofs presented in Section [H estimates on the derivative of the spatial entropy are a 
key element. Since 

and thus 

dt^ T'it) 



^■/^£^-^^(-<->™)^' 



l^T{t) kBT{t) 



it holds 

dE{t) 
dt 



h'{f)dtffiT{t)+ h{f)dtHT(t) 
M Jm 
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h'if) dt^p + ^^ j^ [h{f) - h'if) /] {y - {v)^^,^) fiTit) 

= -ksTit) 1^ h"{f) |V/|^ MT(t) + -^^^ j^ HI) h'if) f] (V (y)^(,)) MT(t). 

In view of Theorem [U the above manipulations (in particular the integrations by parts) can be 
performed, at least for t small enough. The proofs of convergence are based on the following 
strategy: The first term in the above equality is bounded by ~pE{t) for some p > using a 
functional inequality (logarithmic Sobolev inequality (LSI) or a Poincare inequality, depending on 
the functional setting) , while the remainder is shown to be small for 7 small enough thanks to the 
temperature derivative factor. 

2.2.1 Relative entropy estimates 

This case corresponds to the choice h{x) = a; In .t — a:; + 1. In addition of Assumption [2 we suppose 
that 

Assumption 3. There exists an interval I^^^ = [T^f^i, T^J^] (with T^g > and T^g < T* < 
^maxJ such that the family of measures {/^tIts/^^^ satisfies a logarithmic Sobolev inequality (LSI) 
with a uniform constant 1/p, namely 



[ Kf)pT <- f 

Jm P Ja 
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A LSI holds for instance in the following cases: when the potential V satisfies a strict convexity 
condition of the form Iless(y) > vld with v > Q (as a special case of the Bakry and Emery 
criterion [3]), or when the measure /iy is a tensorization of measures satisfying a LSI (see Gross [9]). 
Moreover, when a LSI with constant p is satisfied by Z'^ e"^(«) dq, then Z-;j_^^e-(^(«)+^(«» dq 
(with W bounded) satisfies a LSI with constant p = p c™^ w'-sup w _ 'pjjjg property expresses some 
stability with respect to bounded pertubations (see HoUey and Stroock |12|V In particular, a 
LSI is verified for the canonical measure associated with smooth potentials on a compact state 
space (as is the case here). The uniformity of the constant can be ensured by restricting the 
temperatures to a finite interval isolated from 0. 

Theorem 2. Consider a potential V and an observable A such that Assumptions\^ and\3[hold, and 

an initial data (T°, V'°) with ifj^ e H^(A^), V" > 0, / f/;" = 1, and associated entropy 8{Q) < £* , 

JM 
where 

e* = inf |i(r,A ^ - ^*)^ i(T.i, - ^*)^ \{t]Z ^*)^ i(rLsi _ y*)2 

Then, there exists 70 > such that, for a/Z < 7 < 70, (|15|) has a unique solution (T", t/i) G 
C^([0,r],M) X C*'([0, r], II^(A^)) for all t > 0, and the entropy converges exponentially fast to 
zero: There exists k > (depending on 7J such that 

£{t) <£(0)exp(-Kt). 

In particular, the temperature remains positive at all times: T{t) > min(Tj^fjJ, T^A^) > 0, and it 
converges exponentially fast to T* . 



Remark 2. The proof shows that the convergence rate k is larger when (i) £{Q) is lower (the 
dynamics starts closer from the fixed point and/or closer from a spatial local equilibrium); (ii) a is 
larger (the slope of the function T h^ {A)t is steeper around the target value T* ); (Hi) p is larger 
(the relaxation of the spatial distribution of configurations at a fixed temperature happens faster). 

2.2.2 I? estimates 

This case corresponds to the choice h{x) = -{x— 1)^. In this section, in addition of Assumption [2] 
(and instead of Assumption [3]) , we suppose that 

Assumption 4. There exists an interval /P°in<=arc ^ [yPoincarc^ 1'^™^'=^'^'=] (with TP°;"'=^'^'= > and 
^mhi"''''" < T* < T^™'''^^'') such that the family of measures {/XTlTe/P"'-^"" satisfies a Poincare 
inequality with a uniform constant 1/ p, namely 



M 



h{f)pT<- f |V/|Vt. 
P Jm 



Recah that LSI imply Poincare inequahties (see for instance [5]), so that the second part of 
Assumption m is less demanding than the second part of Assumption [31 and is verified for instance 
for smooth potentials on compact state spaces. A result analogous to Theorem [5] can be stated 
upon defining another bound on the initial entropy (recall that bounds in L? entropy are more 
demanding than bounds in relative entropy [2]). 

Theorem 3. Consider a potential V and an observable A such that Assumptions^^ and^hold, and 
an initial data (r°, V'°) with ■0° G 11^(^1), -0° > 0, / V^ = 1; and associated entropy £{0) < £* , 



JM 
where 

c* • f J [1^^ Tn*\2 (^^ rp^\2 /T-iPoincarc rp^\2 /rriPoincarc T^*\^ 

"1 2 "" '2'' max ^ I '2^ """ ' '2^ max ^ ) 

Then, there exists 70 > such that, for all < j < jo, ()15p has a unique solution {T,tp) G 
C^([0,r],K) X C*'([0, r], II^(A^)) for all t > 0, and the entropy converges exponentially fast to 
zero: There exists k > (depending on 7J such that 

£{t) <£(0)exp(-Kt). 

In particular, the temperature remains positive at all times: T{t) > inm{T^-^f"^'^^^'° , T^^^^) > 0, and 
it converges exponentiallu fast to T* . 

3 Numerical results 

3.1 Single trajectory discretization of the dynamics 



We present in this section a possible discretization of the dynamics (|T3l) . Recall however that the 
general discretization method presented here may be adapted to any dynamics consistent with the 
canonical ensemble. A first obvious strategy is to approximate the expectation by some empirical 
average over K replicas of the system: 

r dq^ = M-^pf dt, 
V 1< fc < A', ' 



dpfe =: ^VViq^) dt - ^M'^Pt dt + ^y2^kBT{t) dW^ 



t ' 



(where {Wf)k=i,...,K are standard independent SA^-dimensional Brownian motions) interacting 
only through the update of their common temperature: 



K 

T'{t)^^^Y.^(^t,Pt)- 
fc=i 
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However, it is often more convenient from a practical perspective to replace the average over many 
replicas simulated in parallel by an average over a single realization. To this end, the dynamics p3|l 
is approximated by the following dynamics: 



' dqt = M^^ptdt, 
dpt = -VV{qt) dt ~ iM-^pt dt + ^2^kBTtdWt. 



dTt^ 



( 



A{qs,Ps)STt^T, ds 



\ 



(16) 



^Tt -T, 



ds 



dt. 



where Wt is still a standard SA^-dimensional Brownian motion, and S is a Dirac mass. The estimate 
of the update term is performed using a trajectorial estimate of the canonical average of A, using 
only the previous configurations sampled at the same temperature. The difference between the 
single trajectory discretization and the many replica implementation is illustrated in Figure [TJ 



T(t):< — ^f^^< >^^ 



X X H — X — )5(X X X — >e< X- 



£ T(t+dt) )i< )0()< — >^ X M^A x|- 



^< — X — X X X X X 



T*- 



Configuration space 



T>T* 



T<T* 




Time 



Figure 1: Left: Discretization using several replicas of the system interacting only through the 
temperature update term. All the replicas are at the same temperature, and estimates of the tem- 
perature update term are obtained from the instantaneous distribution of configurations. Right: 
Single replica discretization. The estimate of the update term at a given temperature is done 
in analogy with the many replica case, using the previous configurations sampled at the same 
temperature. 

The dynamics p6p is discretized using a splitting procedure, with the decomposition 




Ar^ptdt, 



-VV{qt) dt - CM- V dt + ^2^kBTtdWt, 



dqt 
dpt 



0, 
0, 



dTt^ — 



A{qs,Ps) Sxt-T, ds 



/ STt-T, ds 

Jq 



\ 



dt. 



First, the dynamics at a fixed temperature is integrated using some numerical scheme of the general 
form (g"+^,p"+^) = ^ At ,T" {q" , p"") , where (9",^") denotes an approximation of a realization of 
iqnAtiPnAt)i and the current temperature T" is a parameter of the numerical scheme. For instance, 
the BBK scheme [Bl with the modification of [IH] may be used for Langevin dynamics {i = 1, . . . ,N 



11 



indexes the particles): 

n+l/2_ 



^AtJ 



Pi 



n 



-,n+l 



„n+l _ 



At 

At 



-V,.V^(g") - e— + -V2CAtfcBT"Gr 



p. 



n+l/2 



= 1 






2+1/2 



m,. 



At 



- ^V,.V^(q"+i) + -v/2CAtfcBr«Gr 



where {G"}i_„ are independent standard 3-dimensional Gaussian random vectors. Of course, there 
many other implementations of the Langcvin dynamics (see for instance |191 I20j and references 
in [7]). The temperature is then updated after computing the conditional average of the observable 
A along the trajectory. Finally, the algorithm reads 



r ((7"+\p"+i) 



rpn-\-l 



( 



rpn 



m^O 



A(r,P'")xAT(r" 



rpf-{ 



\ 



m=0 



XAT(r" - T") 



7At. 



(17) 



/ 



The ratio in the update term of the temperature is always well-defined provided xat > and 
Xat(O) > 0. The functions xat are approximations of Delta functions, and correspond to some 
binning procedure. For a given temperature, the corresponding bin in the histogram is sought for, 
and the average value of the observable for this temperature as well as the total number of visits 
to this bin are updated. More precisely, discretizing the temperature space into regions of width 
AT, a given temperature T can be rewritten as T = E(T/Ar) AT + Rt^ where E(.t) denotes the 
closest integer to x and — 1/2 < Rt < 1/2. Then, 



Xat(T™ - T") 



1/AT whenE(r™/AT) = E(r"/Ar), 
otherwise. 



The width AT controls the precision of the final result. It should not be too small, so that the 
estimates of the (trajectorial) conditional expectations are performed over enough configurations. 
It may however be refined during the simulation, once close enough to the target temperature. 

The choice of parameters, in particular the parameters of the sampling dynamics (such as Nose- 
Hoover mass or Langevin friction) and 7 required for the update of the temperature, is discussed 
in a specific case in Section 13.2.11 

3.2 A test case: The computation of Hugoniot curves 

We present an application of the general sampling algorithm presented in Section 13.11 to the 
Hugoniot problem described in Section 11.3.11 in the case of noble gases such as Argon. This 
problem was already studied in [15j , where reference results obtained from shock wave simulations 
arc reported. 

The interactions within noble gas atoms arc well-described by a Lennard- Jones potential: 



^(gi,. 



,m) 



\<i<j<.N 



m 



*'^M(7)"-(7) 



In the case of Argon, e/fce — 120 K and a ~ 3.405 A. It is convenient to use dimcnsionlcss 
parameters to present the results, as is done in [T3]. The reduced unit of distance vq = 2^1^ a 
corresponds to the minimum of u, the unit of mass and energy are respectively the mass of an atom 
ra and e. Therefore, the unit of time is r = rt^^Jraje. For Argon, tq = 3.82 A, m = 6.64 x 10^^^ kg. 
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e = 1.66 X 10^^^ J, and t = 2.42 x 10^^^ s. The cut-off used for the computation of the forces is 
i?cut = 2.5(7, and we chose AT = 0.2 K (AT = 1.67 x 10^'^ in reduced units) for the temperature 
histograms. 

The results presented here arc obtained for a crystal reference state, in the cubic face centered 
geometry, at Tq = 10 K with initial density po = 1.806 x 10'^ kg/nr^. The density is chosen so that 
the initial pressure Pq — 0. 

3.2.1 Choice of parameters 

In the case of Hugoniot curves, it is convenient to rewrite the parameter 7 arising in the update 
equation for the temperature (|16p as 

V 

7 



where z^ > is a frequency. The temperature update can then be recast in a form involving only 
dimensionless quantities: 






^(^)=-#^-rf^> A(T)=^o 



5t-Ts ds 


since the observable Ac (or Axx,c)i defined by ([S]) is an extensive quantity homogeneous to an 
energy. This is then also the case for the average quantity At{T). The energy E-^ci ~ k-QT^^i is 
some reference energy per particle. Wc first discuss the choice of the reference temperature, and 
then the choice of the frequency u. 

The other parameters have standard values in molecular dynamics simulations: the time step 
Ai ~ 0.001 in reduced units, while the Langevin friction ^ is such that ^/m ~ 1/t. Our simulations 
for Argon were performed with Ai = 2 x 10~^^ s~^ and £,/m = 10^^ s~^. 

Reference temperature. The reference temperature is computed from the initial configuration 
using the estimator of the temperature derived in |15| . Decomposing the microscopic observables 
associated with the energy and the pressure as 

H = i^kin + ^pot , P = P\i\n + -Ppot , 

with 

the Hugoniot problem Q implies for a given compression c: 
1 2 



T = 



(WldcTo - {V)c\Vo\,T + \{{Pvot)cV„T + {P)\V\o..n){l - c)|I?o|) 



Nk^ 4 - c 
Therefore, an estimator of the temperature in terms of the positions q of the particles only is 

T{q) = ]^ J^ (WiCohTo - V{q) + \{P^ot{q) + {P)\V.ln){l ^ c)|I?o|) , (18) 



which is by construction such that \T\ =T. The reference temperature is chosen to be 

T,ef = f (g°), (19) 

where g*^ denotes the initial configuration. Usually, this initial configuration is obtained by com- 
pressing a perfect lattice in the chosen direction(s). When the compression is not isotropic, the 
estimator ([T5| for the initial temperature should be replaced by a similar estimator where the 
isotropic pressure observable P is replaced by the observable associated with the Pxx component 
of the pressure tensor. The decomposition into kinetic and potential parts of Pxx is done as for P. 
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Typical frequency. Once the reference temperature is set, it is possible to estimate the typical 
magnitude of the normalized observable 



Aciq,p) 



M<1,P) 
NkBT, 



rcf 



for instance by performing some short canonical samplings at fixed compressions and tempera- 
tures. The quantity A^x.c is defined and estimated in a similar manner. In the case of Hugoniot 
sampling, Figurc[2]presents the distribution oi Axx,c for canonical samplings performed at different 
temperatures and anisotropic compressions. The typical values oi Axx,c are of order AAxx,c ^ 0.5. 
This gives a typical size for the term appearing in the equation on the temperature since the tem- 
perature is updated by a factor of order AAxx,c'^At at each time step. The magnitude of v should 
be chosen such that the evolution of the temperature is not too fast (as suggested by Theorems [U 
and [3]). For instance, 

^ AAxx,ci^^t^lO-^ -IQ-'^. (20) 



Tr 



rcf 



Typical time integration steps are At =10 '^ — 10 ^ in reduced units. The parameter i/ can then 
be found from the above heuristic rule. 



0.2 


J 


f-^" 


x" 

X 

< 
V -0.2 


1 


A 








-0.4 




k"^ 






-*- 
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Temperature 



40 



Figure 2: Distribution of Axx,c (in reduced units) for an anisotropic compression c = 0.62 as 
a function of temperature. The error bars are equal to twice the standard deviation of Axx,c 
obtained for a system of iV = 4, 000 atoms, and arc meant to give an order of magnitude of the 
width of the distribution of values of Axx.c- 

In the Hugoniot example, the rule ((20)) gives an estimate of possible values of the frequency. 



For Argon at a compression c = 0.62, 



10"'^^ s ^ seems a reasonable choice. The results of 



Figures [3] and [H obtained for systems of A^ = 4, 000 and N — 32, 000 atoms respectively, show 
that the temperature converges with the method presented here, and that frequencies around 
V = lO""^^ s~^ are indeed a reasonable choice. For smaller frequencies, the convergence is slower, 
whereas larger frequencies trigger fast initial oscillations which may lead to numerical instabilities 
in the scheme. In any case however, the final result does not depend on the value of ly, which 
shows the robustness of the method. It may be interesting to increase the value of v as the typical 
values of the observable decrease during the simulation (see also Remark [Ij . 

3.2.2 Hugoniot curve 

We present in Figure [5] the Hugoniot curve obtained for a system of A^ = 4, 000 particles, using 
the parameters given in Section 13.2.11 with ly = 10"'^^ s~^ for compressions c < 0.7 (temperatures 
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Figure 3: (color online) Plot of the temperature as a function of time (in reduced units) for different 
values of the frequency v (in s^^), for a system of size A^ = 4, 000. 
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Figure 4: (color online) Plot of the temperature as a function of time (in reduced units) for different 
values of the frequency v (in s^^), for a system of size N = 32,000. 

around 4 and less, in reduced units) and v = 10^^ s~^ for higher compressions. This curve is 
obtained by considering many different anisotropic compressions, and computing the temperature 
as well as the associated average pressure P in the system. The results arc in very good agreement 
with [T^. which validates the method. 

4 Proofs of the mathematical results 

4.1 Proof of Theorem [T] 

Existence of a time t > where the temperature remains positive. Assume that V-'(i, ) € 
L^(A^) and ?/; > 0, so that ||V'(^, OIIli = 1- Then, since the temperature derivative is bounded by 
7||A||oo, there exists a time r > such that |T(i) - T°| < r"/2 for all t e [0, r]. For instance, 



yO 



T = 



27PII 



>0, 



(21) 



where ||A||oo = sup{ |^((?)|, (7 £ A^ }. This choice ensures that the temperature remains positive. 
We will show below that il) e L^(A^) since tA € L^(A^) and |A^| < +00. 



Construction of a solution through a fixed point appUcation. We show the existence of 
a solution using a fixed-point theorem. Denote 

ii; = Ci([0,r],R)n| — <T< — ln{T(0)=T"}n{|T'|<7||A|U}. 
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Figure 5: (color online) Hugoniot curve for Argon (in reduced units), computed using the method 
presented in this work (black solid line and crosses). The reference results from [15j are also 
reported (red circles). 



This set is a bounded closed convex subset of C^([0,t],R). Consider the mapping 



9 ■ 



E -^ CO([0,r],H2(X)) 



C^{[0,T],R)r\E 



E 
9{T) 



T ^ Vt ^ 9{T) 

The construction of V't from a given temperature function T is performed using the equation 



(22) 



dti^T = kBT{t) V 



fJ'T(t) 



\f^T{ 



= fcBr(i)AV'T + vy • VV'T + A^^V'T, (23) 



and the temperature function g{T) is then recovered using the second equation in (jlSp . namely 



Jo Jm 



dqds. 



Regularity of ■07- given T. The equation (P5|) is of the general type 

dtilT + CT{t)lpT = 0, 

where the family of operators {CT{t)}t>o have constant domains D ~ H^(A^), and 



(24) 



CT{t)i^^-kBT{t)W 



f^Tit)"^ 



i, 



f^T(t) 



= -V • (fcBT(i) vtA + vy^) 



It is enough to show the regularity of the solution of dtip +{''] + CT{t))ip = since the latter 
equation is related to the original equation ((M|) through the transform ipx = exp{rit)ilj. 
For a fixed t and / S L^(A^), consider the equation 

iv + CTit))u = f. (25) 

Lax-Milgram's theorem shows that there is a unique solution u e H^(A^). Indeed, {u,v) i— > 
{u, (r] + CT{t))v) can be extended to a functional on H^(A^)^ using 



\{u,{^ + CT{t))v)\ 



/ j-juv + kBT{t)Vu ■ Vv + vVu ■ W 

JM 



<C\\uU.\\vU., 
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where the constant C depends only on ||VF||oo and on T° (actually on the upper bound of the 
temperature). Besides, u £ D i-^ (m, (77 + CT{t))u), extended to a quadratic form on H^(A^), is 
coercive on this space for i] large enough since 



{u,{r] + CTit))u) = / rju^ + kBT{t)\Vu\^ + uVu-VV 
Jm 

> M!! / |V^r + W u'-^\\VV\\^[a I \Vu\-' + - 



M 



M 



M 



a 



M 



> c\\u\\ii^, 



for r] large enough (choosing first a < kBT^/W^VWoo and then 77 > ||Vyj|oo/2a). Again, the 
constant c depends only on ||Vy||oo and on T°. The solution u of ([25|) is in fact such that 



Am = 



kBT{t) 



{r]u -VV-Vu- uAV - /) E L'^{M), 



so that u E D. It is easily seen that the mapping f 1-^ u = {X + t] + Crit))^^ f for A G C with 
Re(A) > is continous, and bounded by M/(l + |A|) for some M > 0. 
Besides, for < r, s, t < r, 

II [(r? + CTir)) -iv + Cris))] (v + CTit)r'\\ciL^,L^) = fcB|T(r) - T(s)|||A(r; + CTit))-'\\ciL^,Ln 

< c\Tir)~Tis)\<C\r~sl 

for some constants C, c > since T is C^. 

The results of Kato [Mj then show the existence of a unique solution t/jt E C°([0,t], H^(A^)) 
and hence ipr G C'dO, r],H^(A^)). Bounds on the H^(A^) norm oi ipr in terms of the bounds on 
T are obtained with a Gronwall inequality: 



dt [ ^WMlw 



[iPt - AiPt)V {kBT{t)V'ipT + VVipr) 



M 

-kBT{t) f IWtP + IAtAtI' 

JM 



AlpT^V ■ VljJT 



M 



Jm Jm 



Using the following Cauchy-Schwarz inequalities to bound the terms involving AipT' 

fee J 



A^JJT^V ■ \/ijT 



M 



< M! / |A^,|. . IIV^II^ 



M 



kBT° JM 



iVVrr, 



iPtAiPtAV 



M 



< 



fcRT" 



lA^Tp + 



M 



l|AT/||g. 

fcRT" 



\^i 



M 



it follows 



dt(l:UT\\H^] < -fcBr(t) / |V^Tp+ / I^TpATZ + ^rVV^-VV'T 
V^ / Jm Jm 



\yv\\l + \\Av\\ 

kBTO 



IV'tII^i 



where the constant K > depends only on V and its derivatives, and on r° (actually, the lower 
bound of the temperature) . This shows that 



||7ZT(t,-)||Hi<exp(ri^)||V'"||Hi 



for all < t < T. 
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Regularity of g{T) given ij}. From the equation on the temperature: 

(7(r)'(t) = -7 / A{q)i'T(t,q)dq, 
JM 

it foUows 



g{T)"{t) = 7 / {WA-WV- kjiT{t)AA)i,T, 

JM 



and 



g(T)(3)(i) = -7 /" V • (V^ • VF - kBT{t)AA){'K7V^/jT + kBT{t)\7i'T) - iknT'{t) f AAi/'t- 

JM JM 

Therefore, \giT)'^^\t)\ < ^CAy\\T\\c^il + \\T\\ci)\\Mtr)\\m, ^i^ove Ca^v depends only on \\M\ 
and on L°° bounds on A, V and their derivatives (up to the order 3 for A, and to the order 2 for 
V). Uniform bounds on g{T)" are recovered by time integration: 

sup |g(r)"W| < |g(r)"(0)|+r||.g(T)(-^)|U 

Q<t<T 

< -f{CA,v\\^°\\H^+TCA.v\\T\\c^il + \\T\y)\ma^i[o.T]Mn) 

< 7pAy + TCA.y||r||ci(l + ||T||ci)cxp(AV)] |K||hi < M, 

for some M > (depending on r, ?/>o, T°). Since ||.g(T)'||oo < 7PI|oo, bounds on g{T),g{T)' 
consistent with the set E are still valid with the choice of r given by (PT|) . This shows finally that 
g[T) belongs to a bounded set of C2([0,t],R). 

Existence of a fixed point. The mapping g is continuous from the bounded closed convex 
set E to istself, and in fact compact since the injection C^([0,r],IR) n {\T"\ < M} (1 E ^^ E is 
compact. Schauder's theorem (see for instance [2T]) shows that there exists T G E such that 
g{T) ~ T, therefore giving the existence of a solution to the equation (fT^ . 

Uniqueness. Consider two solutions (Tijipi) and {T2,ip2) starting from the same initial condi- 
tion (r°, -0°). Then, on the interval [0,t], 

4('im(i)-T2(t)f) = -7(Ti(i)-T2(t))^A(Vi-V2) 



dt \2 



< ^\M\\\A\\^ |Ti(t) - T2{t)\ 1101 - i'2\\L-. 



Moreover, 



_^^TU0^^2W f v(0i + V2) • V(01 - ^2) 

^ JM 



< Ci||0i-02||£2+C2|Ti(i)-T2(t)|M |V(Vl+02)r, 

JM 

for some constants ci, C2 > 0. Indeed, the the V('0i — 1P2) contributions of the first and third terms 
on the right hand side of the first line can be cancelled using the inequality 2ab < 0^/1] + ■qlP' and 

_,^ Ti(t)+T2(t) r |v(0,_^,)p<_^/ |v(^i-^2)r 



M 2 j^ 
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In view of the uniform bounds on ||V'i/'i(i, OIIl^, ||VV'2(ij OIIl^, 
for some C > 0. As 

^(ri(o) - T2(o))2 + ^11^-1(0, •) - v^2(o, OIIl^ = 0, 

a Gronwall inequality shows that Ti = T2 and V-'i = "02 for all < t < r. 

Positivity of ijj for relative entropy estimates. It is important that -0 > for the Fisher 
information to be defined in the case of relative entropy estimates (see the proof of Theorem [5] in 
Section [4.2p . In any cases, the density "0 has to remain non-negative. 
Consider 

^t,q)=i^{s{t),q), (26) 

where the function t ^—t s(t) is defined through the ordinary differential equation: 

1 



s'it) = 



kBT{s{t)y 



The function s is well defined (since T > 0), continuous and in fact C^, and strictly increasing 
hence invcrtible. The evolution of the density $ is obtained by solving 

3,^ - A<1. + ^-^V . ($VT.) = A* + ^V . i^VV), (27) 

with T{t) ~ T{s{t)). The density '0 is recovered with (|26p since s is invertible. The function 
= *&/,/A'f(t) satisfies 

dtcj){t, x) = A0(t, x) + W{t, x)(t){x), 

where the potential 

W{t,x) = — L— I Ay(x) ^|v\/(x)p - Z^V{x)\ 

2knTit) \ 2fcBT(i) T{t) 'J 

is bounded. Since fix is positive and uniformly bounded from below, it is enough to show that 
4>{t, •) > to have ip{t, •) > 0. The function <f> = cxp(f||IF||oo)0 satisfies 

[dt - A)0 = iWit,x) + ||I^||oo)0 > 0. 

Therefore, 0(t, •) > inf0(O, •), and so, 4>{t,x) > mexp(— t||M^||oo) > provided 0(0,-) > to > 0. 
This shows that -0 > provided 0" > 0. Similarly, > provided V'° > 0. 

4.2 Proof of Theorem [2] 

General strategy of the proof. We first consider the solution of (fT5|) on a time interval [0, r] 
as given by Theorem [I] possibly reducing r so that T{t) e /^ n T^^^ for all < t < t. We then 
show (see below) that this implies £'{t) < ~K£{t) for some some k > 0. This will show that £ 
decreases, and so T remains in I^ f] T^si since \T{t) - T*\ < ^2£{t) < ^2£(0) < V2£*. The 
solution can then be continued at all times t > using repeatedly the arguments of the proof 
of Theorem [1] upon replacing the set E in the definition of g given in ([22]) by the closed convex 
bounded set E = G\[nT, (n+l)T],R)n{r G I;^nT^^^}n{T{nT) = r"^}n{|r'| < 7PII00}, where 
T"'^ is the final value of the solution on the interval [{n — 1)t, tit] (for 71 > 1). The main difference 
with the proof of Section 14.11 is that uniform (upper and lower) bounds on the temperature are 
provided thanks to the decrease of the entropy. 
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Proof of the Gronwall inequality £'{t) < —K£{t). With the function h chosen here, 



dE{t) 
dt 



-knT{t) 



|V/|^ 



T'{t) 



^ / ^^(*) k^ntyjM 



u-^)[y-{y)nt))mt) 



< -pkBT{t)E{t) + 

< -pkBT{t)E{t) + 



2\T'{t)\\\V\\ 
ksTity 

MT'{t)\\\v\\ 

kBT{ty 



I l/-l|MT(t 
J M 



'-Vm. 



where we have used the Csizar-KuUback inequahty 



M 



|/-l|MT(t)<2/B(i) 



in the last step. The derivative of the temperature can be bounded using 

T'{t) = -iJ^Mt = -7 ((^>T(t) + /^ ^(/ - l)/^^(*)) ' 



which shows that 



\T'it)\ < 7(|(^)T(t)-(^)i 

< l(^a\T{t)-T*\ + 

< ^(^a\T{t)-T*\+2\\A\U^/E{t) 



Aipt - / ApT{t) 

M J M 



Atpt - / AfiT{t) 

M J M 



using Assuniption[2]and the Csizar-KuUback inequality to bound the second term on the right-hand 
side. These estimates are also helpful for bounding the derivative of the temperature entropy: 



d 
dt 



*\2 



'-{T{t)-T*) 



< -7 ((^)t(0 - (^)t*) {T{t) T*) + 27||A||^ \T{t) T*\ y^W) 

< -a-/{T{t) - T*f + 27||A||oo \T{t) - T*\ ^fW). 



where Assumption [2] has been used again. Gathering the estimates, 



d£{t) 
dt 



< 



< 



< 



-pkBT{t) E{t) 



^l\\V\\o 

kBT{tf 



-VW){a \T{t) - T*\ + 2\\A\\^^E{t)) 



aj{T{t) - T*f + 27PII00 \T{t) - T*\ ,/Eit) 
pkBTit) ^l&Mk) E{t) aj(T{t) - T*f 
4a 11^11 



7(2||A||„, 

pkBT{t) 



a-rj{\\A\ 



kBT{tY 

87rHoo||^ll. 

kBT{t)^ 

2a||yi| 



y/W)\T{t)-T* 
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\A\\ 



kBTity 
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where the inequality ab < — a -\ — b (for any r/ > 0) has been used in the last step. This estimate 

2ri 2 

shows that, with j] — ^/7 for instance, and choosing 7 small enough such that 
there exists k > such that 

^ < -<')■ 

Therefore, the entropy decreases, and the initial bounds on the temperature are satisfied at all 
times. Gronwall's lemma finally shows that £{t) -^ exponentially fast. 

Remark 3 (Nonlinear feedback). Convergence results for a nonlinear temperature feedback of the 
form p4p for a nonlinear function f such that f{x) = if and only if x = 0, can be proved by 
extending the proofs of this section to the case when 

f(x) 
Vi? > 0, 3ci, C2 > such that Vx e \-R, R], a < ^^-^^ < C2. (28) 

X 



The constant R is related to the maximal value of / Aipt , which is dictated by the initial entropy 

JM 
£(0).' if £{t) < £{0), then the temperature T{t) belongs to a bounded set Ism), and so, for usual 

entropy functions, 

f A{q)MQ)dq <\{A)T(t)\+ I A\f-l\lJiT(t)dq< sup \{A)t\ + c^£{t) < R. 
JM JM TG/£(o) 

where R depends only on £{Qi), and c depends on the choice of the entropy function (c = 2 for 
relative entropies). This shows why the bounds (|28p on some compact interval only are sufficient. 
For the choice f{x) ~ x + x^ for instance, ci = 1 and C2 = 1 + R^ . Such choices may be helpful 
in numerical simulations to accelerate the convergence toward the temperature of interest at the 
early stages of the process. 

4.3 Proof of Theorem [1 

The proof follows the same lines as the proof of Theorem ^ and we present only the modifications 
required in the case considered here. The first change in the proof is in the bound of the derivative 
of the spatial entropy E{t). It holds 

^ = -^-B^w L i^^i^ ^-(*) + iBw L t'^^^ - ^^ - '^ ^^ (^ - ^'^^-<*^) ™ 

= -pk^T{t) E{t) - -^^ j [h{f) + f-i]{y- (F)^(,)) ^iTit) 
. ( , r(A '^\T'{m\v\\oo \^u. , 2|r(t)|||y||oo [ ,, ^1 

2|T'(i)|||y||oo\ ^,,, , 2V2\T'm\V\\ 



. -W^^m'^-^^^^^^)m. -^l^'^;^^- Vm 



where we have used the Cauchy-Schwarz inequality 



\f-i\mt)<\ IZ-iIVtwa// tiTit) = vmt)- 

M V J-M V J-M 
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A Cauchy-Schwarz inequality is also used to establish a slighty different bound on the temperature 
derivative: 

\r{t)\ < 7 (a \T{t) -T*\ + \\A\\^y^2m) , 



and 



d 
di 



*\2 



(r(i) -T*) 



< -aj{Tit) - T*f + 7||A||oo \T{t) - T*\ yj2E{t). 



The spatial entropy derivative can then be bounded as 
dE{t) 



dt 



< 



pkBT{t) - 7 



M\A\loo\\V\\ 

fcBT(t)2 



f 



E{t) 



E(')+ ^rJ^^!f \Tit)-n {E(t) + vmi 



For the total entropy, 
d£{t) 



dt 



< - 



pkBT{t) - 7 



+V2- 



7 



2a\\V\\ 



M\Moo\\v\\c 

kBT{tf 



E(t) 



E{t) - a-f{T{t) - T* 



E{t) 



\A\\ 



\T{t)^T*\./W)- 



kBT{tY \" V 2 

This estimate is analogous to the one obtained in the relative entropy case, except for some 
additional factors ^J E{t), which are bounded by yj£{t), hence by \/£{0) when the total entropy 
decreases. However, it is easily seen that, taking initial conditions £(0) < £*and for 7 small 
enough, the function £ is decreasing and bounded by £{Q). There exists then k > such that 
£' [t) < —K£{t), which again implies the longtime convergence. 
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